#include <math.h>
#include <stdint.h>
#include <math-barriers.h>
#include <math-narrow-eval.h>
#include <math-svid-compat.h>
#include <libm-alias-finite.h>
#include <libm-alias-double.h>
#include "math_config.h"

#define N (1 << EXP_TABLE_BITS)
#define InvLn2N __exp_data.invln2N
#define NegLn2hiN __exp_data.negln2hiN
#define NegLn2loN __exp_data.negln2loN
#define Shift __exp_data.shift
#define T __exp_data.tab
#define C2 __exp_data.poly[5 - EXP_POLY_ORDER]
#define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
#define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
#define C5 __exp_data.poly[8 - EXP_POLY_ORDER]

/* Check if the result overflowed to infinity.  */
static inline double
check_oflow(double x)
{
  return WANT_ERRNO ? __math_check_oflow(x) : x;
}

/* Check if the result underflowed to 0.  */
static inline double
check_uflow(double x)
{
  return WANT_ERRNO ? __math_check_uflow(x) : x;
}

static inline uint64_t
asuint64(double f)
{
  union
  {
    double f;
    uint64_t i;
  } u = {f};
  return u.i;
}

static inline double
asdouble(uint64_t i)
{
  union
  {
    uint64_t i;
    double f;
  } u = {i};
  return u.f;
}

static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki)
{
  double_t scale, y;
  if ((ki & 0x80000000) == 0)
  {
    sbits -= 1009ull << 52;
    scale = asdouble(sbits);
    y = 0x1p1009 * (scale + scale * tmp);
    return check_oflow(y);
  }

  sbits += 1022ull << 52;
  scale = asdouble(sbits);
  y = scale + scale * tmp;
  if (y < 1.0)
  {
    double_t hi, lo;
    lo = scale - y + scale * tmp;
    hi = 1.0 + y;
    lo = 1.0 - hi + y + lo;
    y = (hi + lo) - 1.0;

    if (1 && y == 0.0)
      y = 0.0;

    if (sizeof(
            ({
              __typeof(0x1p-1022) __x;
              if (sizeof(0x1p-1022) <= sizeof(double) || __builtin_types_compatible_p(__typeof(0x1p-1022), _Float128))
                __asm(""
                      : "=x"(__x)
                      : "0"(0x1p-1022));
              else
                __asm(""
                      : "=t"(__x)
                      : "0"(0x1p-1022));
              __x;
            }) *
            0x1p-1022) <= sizeof(double) ||
        __builtin_types_compatible_p(__typeof(({
                                                __typeof(0x1p-1022) __x;

                                                if (sizeof(0x1p-1022) <= sizeof(double) || __builtin_types_compatible_p(__typeof(0x1p-1022), _Float128))
                                                  __asm(""
                                                        : "=x"(__x)
                                                        : "0"(0x1p-1022));
                                                else
                                                  __asm(""
                                                        : "=t"(__x)
                                                        : "0"(0x1p-1022));
                                                __x;
                                              }) *
                                              0x1p-1022),
                                     _Float128))
      __asm __volatile(""
                       :
                       : "x"((
                                 {
                                   __typeof(0x1p-1022) __x;
                                   if (sizeof(0x1p-1022) <= sizeof(double) || __builtin_types_compatible_p(__typeof(0x1p-1022), _Float128))
                                     __asm(""
                                           : "=x"(__x)
                                           : "0"(0x1p-1022));
                                   else
                                     __asm(""
                                           : "=t"(__x)
                                           : "0"(0x1p-1022));
                                   __x;
                                 }) *
                             0x1p-1022));
    else
      __asm __volatile(
          ""
          :
          : "f"(({
                  __typeof(0x1p-1022) __x;
                  if (sizeof(0x1p-1022) <= sizeof(double) || __builtin_types_compatible_p(__typeof(0x1p-1022), _Float128))
                    __asm(""
                          : "=x"(__x)
                          : "0"(0x1p-1022));
                  else
                    __asm(""
                          : "=t"(__x)
                          : "0"(0x1p-1022));
                  __x;
                }) *
                0x1p-1022));
  }

  y = 0x1p-1022 * y;
  return check_uflow(y);
}

static inline uint32_t top12(double x)
{
  return asuint64(x) >> 52;
}

double __ieee754_exp_sse2(double x)
{
  uint32_t abstop;
  uint64_t ki, idx, top, sbits;
  double_t kd, z, r, r2, scale, tail, tmp;

  abstop = top12(x) & 0x7ff;
  if (abstop - top12(0x1p-54) >= top12(512.0) - top12(0x1p-54))
  {
    if (abstop - top12(0x1p-54) >= 0x80000000)
      return 1 ? 1.0 + x : 1.0;

    if (abstop >= top12(1024.0))
    {
      if (asuint64(x) == asuint64(-(__builtin_inff())))
        return 0.0;
      if (abstop >= top12((__builtin_inff())))
        return 1.0 + x;
      if (asuint64(x) >> 63)
        return __math_uflow(0);
      else
        return __math_oflow(0);
    }
    abstop = 0;
  }

  z = __exp_data.invln2N * x;
  kd = (z + __exp_data.shift);
  ki = asuint64(kd);
  kd -= __exp_data.shift;

  r = x + kd * __exp_data.negln2hiN + kd * __exp_data.negln2loN;
  idx = 2 * (ki % (1 << 7));
  top = ki << (52 - 7);
  tail = asdouble(__exp_data.tab[idx]);
  sbits = __exp_data.tab[idx + 1] + top;
  r2 = r * r;

  tmp = tail + r + r2 * (__exp_data.poly[5 - 5] + r * __exp_data.poly[6 - 5]) + r2 * r2 * (__exp_data.poly[7 - 5] + r * __exp_data.poly[8 - 5]);

  if (abstop == 0)
    return specialcase(tmp, sbits, ki);
  scale = asdouble(sbits);
  return scale + scale * tmp;
}
